Origin of Anomalous Diffusion and Non-Gaussian Effects for Hard Spheres: Analysis 

of Three-Time Correlations 
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We present new simulation results on a hard sphere system at high densities. Using three-time 
correlations, we can account for the anomalous diffusion, which results from a homogeneous back- 
dragging effect. Furthermore, we calculate the non-gaussian parameter and connect it to the existence 
of dynamic heterogeneities. 



I. INTRODUCTION 

In recent years the nature of the non-exponential 
relaxation of glass-forming liquids close to Tg has been 
elucidated by a variety of experimental techniques; 
see and H for an overview. One way of getting 
additional information is to analyse three-time besides 
the standard two-time correlation functions, as realized 
by multidimensional NMR experiments. On this basis it 
is possible to decide whether the non-exponentiality is 
mainly due to a homogeneous scenario (intrinsic non- 
exponentiality related to systematic back- and forth 
dynamics) or a heterogeneous scenario (superposition 
of different exponential processes). In a previous 
publication these concepts have been applied to computer 
simulations of a hard sphere system which experimentally 
is well represented by colloidal suspensions . There we 
observed that the dynamics at short times (/3-regime) 
is dominated by the presence of a cage formed by 
surrounding particles. This results in a systematic back- 
dragging effect, corresponding to a mainly homogeneous 
scenario. In the a-regime most particles have escaped 
their initial cages. Then they encounter fast or slow 
regions in the liquid, arising, e.g., from inhomogeneitities 
in density. The simplest model for this situation would 
be an ensemble of independently diffusing particles with 
different mobilities. 

The goal of the 

present work is to quantify the hard sphere dynamics in 
terms of homogeneous and heterogeneous contributions 
for all time scales. In extension of our previous work, 
we introduce appropriate measures for homogeneous and 
heterogeneous contributions. They help to elucidate the 
physics of the sublinear diffusion in the /3-regime and the 
non-gaussian parameter. 



II. SIMULATION DETAILS 

We perform conventional Monte Carlo dynamics with 
periodic boundary conditions on a system of a thousand 
hard spheres with a size polydispersity of 10 percent. The 
algorithm is described in Mm- For the P and a regimes. 



the mean step size has been chosen to give an acceptance 
rate of 50 percent for each move. We checked, that, apart 
from an overall scaling, the dynamics is not influenced by 
changing the mean step size for these time regions. 

At very short times, we perform separate simulation 
runs with a reduced step size, in order to obtain the 
microscopic dynamics. This is necessary, because we want 
to define a density-independent time scale by fixing the 
short time diffusion constant, i.e. {r'^{t)) = QDot, where t 
is very small. We use Dq = -^ and set the unit of length 
equal to the mean particle diameter, as is done by Fuchs 
et al. II . 

The highest density, analysed in this paper, is = 
(|7ri?f)^ = 0.58. At this high packing fraction, we 
have chosen a very long equilibration time (5 • 10^ 
MCS) to avoid aging effects during the measurement. 
Interestingly, the non-gaussian parameter (NGP) turns 
out to be a very sensitive indicator for the degree of 
equilibration. 




Figure 1 

(a) The two-time correlation function _F2(t) = {cos q{r{t) — f{0))) 
for q = 2iT which is close to the maximum of the structure factor. 
The curves are plotted for the densities (p = 50, 53, 56, and 58%, 
from left to right, (b) The mean square displacements (r^(t)) for 
the densities in (a), again from left to right. 

In figure || we plot the incoherent scattering function 
F2{t) = (cos(f(r(i) — r{0))) for g = 27r and the mean 
square displacement, clarifying the time scales of our 
simulations. 



III. THREE-TIME CORRELATIONS: THE 
HOMOGENEOUS PART 

As shown in our recent work H, we can gain 
new insight into the dynamics by analyzing three-time 
correlation functions. For a given time t we define 
the conditional probabilities 7*11(^12 • ?\ii|?'oi;i), Pxi^u ■ 
uoi\roi;t) and Pzzizi2\zoi;t) with rv„„ = r{nt) - 
f(mt), r„in = I^Vrml- Hcrc uqi is an arbitrary vector 
perpendicular to roi, the hat denoting a unit vector. We 
abbreviate a;i2 = ^12 • ^31 and 2/12 = ri2 • moi for later 
purposes. Thus, p||(a;i2|roi; t) denotes the probability, 
that a particle moves during the second time interval 
the distance a;i2, projected onto the vector foi, under 
the condition that it made a distance roi in the first 
time interval. The probability function p±{yi2\roiTt) 
is defined analogously, yielding information about the 
distance 2/12 travelled perpendicular to the first step. 
Finally, Pzzizi2\zoi;t) includes the information about 
the projection of both steps onto a randomly chosen 
direction. 
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Figure 2 

The conditional probabilities Pzz{zi2\zoi;t), p||(a;i2|r'oi; i), and 
P±{yi2\'>'oi',t) for a time t fs 400 in the /3-reginie, where a;i2 = 
^I2''0l and j/12 = f'i2Uoi. 201 and zi2 denote the projections of 
the subsequent displacements onto an arbitrary direction. The dark 
areas correspond to high probabilities. The average values c^z (201)1 
C||(roi), and c±{roi) are indicated by the light lines. Note that the 
unit of length is greater than in lil by a factor of 2. 



For non-exponential relaxation, as displayed by 
supercooled liquids, memory effects are present, resulting 
in a strong dependence oi Pzz and p\\^± on zqi and rgi, 
respectively. We see the typical situation in figure HI, 
where all these probability distributions are shown for 
a density oi (p ~ 58% at a time t « 400 in the /3- 
region. The main features of these plots are the non- 
zero displacements, Czz{zoi) = (Z12), C||(roi) = (2;i2) and 
c±iroi) = (2/12) as indicated, and the widths azz{zoi), 
cr[|(roi) and a±{roi), where, e.g.. 



O'wiroi) ^ / da;i2P||(a;i2|roi)(a;i2 - C|i(roi)) . 

The possible dependence on t of these quantities should 
be kept in mind. 



As expected from the definition of p±, the value of 
c_l(?'oi) must be zero. We discussed in 0] that in the 
heterogeneous limit, there is no back-dragging effect, i.e. 
Czz{zoi) — C||(roi) = whereas in the homogeneous case 
we have <Jzz.\\. ±{1^01) — const and c^^ ||(roi) = const. 

As shown in [Q, Czz{zoi) - which, by the way, is equal 
to C||(roi) - can account for anomalous diffusion in the 
/3-domain. To be precise, we made use of the fact that 
Czz{zqi) « — czoi, where c = c{t) is a constant. If 
we additionally assume that Pzz{'\zoi) is gaussian with 
constant width azzizoi) = Cj straightforward integration 
yields a connection between c{t) and the approximate 
logarithmic slope 
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of the mean square displacement, namely, 
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This estimate becomes wrong for longer times, because 
the back-dragging is not linear anymore. Mainly this 
comes from the fact that particles start to leave their 
cages, resulting in a more complex behaviour. This may 
explain the deformation of Czz{zoi) and C||(roi), which 
now look more like the first part of a negative sine. To 
account for this effect, we define an 'effective' slope of 
Czz{zqi), namely. 



u\ - J (izoipizoi)zoiCzz{zoi] 
CeS[t) = 



(^0^) 
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where all quantities depend on t. 

Using Ceff(i) instead of c(i) as input for Eq. (1), 
the approximation of f3 becomes exact. With {r'^{t)) = 
3(z2j) = i{zl^), we find 

(^02) = ((^01 + 212)') = 2{{zl^) + (Z01Z12)) 
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If we had the purely heterogeneous scenario, where the 
conditional probability for the second move does only 
depend on the length of the first one, i.e. Pzz{zi2\zoi] t) = 
Pzz{zi2\\zoi\;t), then 

{zoi)ccsit) = / dzoidzi2Pzzizi2\\zoi\)pizQi)zoiZi2 = 0, 

because p(zoi) is an even function. Therefore, Cgfj(t) can 
be regarded as an appropriate measure for homogeneous 
contributions, and in the linear case, it reduces to the 
slope c{t). Hence homogeneous contributions strongly 
influence the nature of the anomalous diffusion. 



IV. THE NON-GAUSSIAN PARAMETER 

A way to approach the high-density features of the self 
part of the van Hove correlation function is to analyze 
the non-gaussian parameter 



a2{t) 
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It measures the deviation from gaussian behaviour and 
consequently has to vanish for i ^ 0, because the 
dynamics on the very microscopic scale is brownian by 
definition. 

The NGP has been calculated for other systems, e.g. 
Lennard-Jones liquids [|| or soft disks Q. In every case, 
a maximum has been found in the /3 region, which is the 
more pronounced, the lower the temperature is. 



et al. have found the same dependence for their binary 
Lennard-Jones mixture 101 . 



In figure IV (b), we see 0^2, poly (0 ^^^ different volume 
fractions in a double logarithmic plot. For small times, 
one perceives a linear ascend which corresponds to 
an exponent of approximately 0.3. Later on, a2.poiy 
reaches its maximum value at a time ta^ in the 
late /3-region and then slowly decreases again. As 
we can see, the maximum value of 02, poly grows 
with density, i.e. 02, max — 0.17,0.24,0.56,2.33, for 
(p = 0.50, 0.53, 0.56, 0.58, respectively. This strong 
dependence indicates a tremendous change of the 
dynamics for densities (j) > 0.56. 



THREE-TIME CORRELATIONS: 
HETEROGENEOUS PART 
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Figure 3 

(a) The non-gaussian parameter 02 = (3(''*})/(5(r^)^) — 1 for 
density cji = 56%, calculated in different ways. For Qj small' '^2 big 
and 02 the average is over the 15% smallest, the 15% biggest 



and all particles, respectively. In a. 



2, poly 



the different sizes are 



treated separately, averaging afterwards, as discussed in the text, 
(b) The NGP a.^ p^jy for the densities </> = 50, 53, 56, and 58%, 
from bottom to top. Notice that the axes in (a) and (b) have been 
chosen logarithmically. 



In figure |rv| (a) we see a^ it) for <^ = 56% as a dashed 
line. It clearly does not decay to zero, but assumes a 
value 02(00) « 0.15. However, we must be careful with 
the definition of a2 , because the polydispersity can cause 
a trivial non-gaussianity, due to a size-dependent particle 
mobility. Therefore, we should calculate the NGP for 
every particle diameter separately, and take the mean 
value over the size distribution afterwards, i.e. 

a2,poly(0 = {a2,R{t))R. 



Figure IV (a) shows this quantity for comparison. 
Thus, at this density, 02, poly (00) = 0. The other two 
curves show the NGP, averaged over the 15% largest 
and smallest particles, respectively. Differing in their 
maximum values by about 10%, they exhibit a time 
separation of their maxima by a factor of three. Kob 



It has been suggested that the value of 02 (t) is 
intimately connected with the existence of dynamic 
heterogeneities. For example. Hurley and Harrowell used 
a model of fluctuating mobilities to account for the non- 
gaussian effects in a two-dimensional liquid P|. However, 
there are a few input parameters for this model, which 
must be adjusted to explain the simulation results, e.g. 
the functional form of the mobility autocorrelation. So, 
it is not clear to what extent the concept of fluctuating 
mobilities meets reality and a deeper understanding 
is still necessary to clarify the underlying physical 
mechanisms. 

Again we analyse three-time correlations, now 
employing the information content of the widths <Jzz.\\.±- 
The important observation in figure HI is that the length 



of the first step of a tagged particle has an influence on 
the mean size of its subsequent step, i.e. cTzzizoi) and 
(T||^(roi) grow with zqi and roi, respectively. This can 
be understood in the following way. A fast particle in 
the flrst time interval on average remains fast during the 
second one, and so, by looking at Pzzi-\zoi) or P||,j_(.|r'oi) 
for large zqi or rgi, we select the most mobile particles. If 
we now calculate the degree, to which azz and cry j^ grow 
with zqi and rgi, respectively, we have a direct measure 
to which extent the dynamics is ruled by heterogeneties. 
For this purpose we define the quantity 



ait) 



(a2}2 



(3) 



where the brackets denote an average over the first step, 
i.e. (..} = /droi..p(roi) or (..} ee /dzoi..p(zoi). Here, a,a 
and p have to be decorated with zz, || or _L, depending 
on the probability function that is being analyzed. (Note 
that a constant width g results in a vanishing a(t).) 
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Figure 4 

The quantities au{t) and aj_{t), as defined in the text, for the 
density <f> = 58%. The diff'erence in the maximum values illustrates 
the fact that heterogeneities at intermediate times possess a 
directionality. For times greater than Tq fa 2 • 10*, both curves 
meet each other, indicating the loss of anisotropy. 

Figure ^ shows a||(t) and a^{t) for the density (f) ~ 
0.58. We see that the broadening with rgi, along the 
direction of the first step, is greater than perpendicular 
to it, resulting in a larger value of a||(t). For longer 
times t > Ta, however, both curves coincide, indicating 
that there is no motional anisotropy with respect to 
the previous step anymore. Collective flow patterns, as 
observed for example by Donati et al. Q, could cause 
such a directionality at intermediate times. 




Figure 5 

The averaged a^^ = (an + 2a^)/3, compared with the NGP 02 
and the effective slope c„g(t) for density </> = 58%. a-^d ^nd 02 show 
the same behaviour with time, their maxima falling into the late 
/3-region. Cgf|:(t), in contrast, is large only for times t < Ta- Note 
that the units are arbitrary. 

It is clear, that a(t) ^ should imply a nonvanishing 
NGP a2{t). Fig. M shows the relation between a2 and 
the averaged a^dit) = ^ — —■ This surprising similarity, 
concerning the time dependence of both curves, suggests 
a deep physical connection. For a molecular liquid, this 



fact has already been reported by Qian et al. jTl) , using 
another quantity similar to a^dit)- Ceff(i) is shown for 
a comparison of dynamical time regimes, reflecting the 
homogeneous contributions to the complex relaxation 
process. 

VI. DISCUSSION 

As a major result of this work we were able to quantify 
the homogeneous and the heterogeneous contributions 
to the non-exponential relaxation. We found that the 
anomalous diffusion is mainly related to homogeneous 
contributions, explicitly expressed by Eqs. (1) and (2), 
whereas the non-gaussian effects, displaying a maximum 
for the crossover between /3- and a-relaxation, are mainly 
related to the heterogeneous contributions. This directly 
shows that for the system, studied in this work, the 
non-gaussian effects, if at all, are only mildly related 
to jump contributions. The maximum value of a2 is of 
the same order as for the case of Lennard- Jones systems 
B. It remains the interesting question whether also for 
the latter a2 is mainly determined by heterogeneous 
contributions. 

Relating the anomalous diffusion to homogeneous 
contributions, as discussed above, is only valid for 
stationary processes. As outlined, e.g., in Ref. [ p^ , 
continuous time random walks, which by definition are 
purely heterogeneous in nature, arc also able to generate 
anomalous diffusion behaviour. This, however, explicitly 
requires non-stationary conditions. 

We would like to mention that a(t), defined in Eq.(3), 
only takes into account heterogeneous contributions, 
related to the roi-dependence of ctjm (^oi) and neglects 
contributions resulting from the tqi -dependence of 
C||(roi). A simple model system with <J\\,±{rQi) = const 
and cii (^oi) oc — roi is an ensemble of diffusive particles in 
identical harmonic potentials . The dynamics is purely 
gaussian, i.e. a2 = 0. Here the heterogeneity is related to 
the fact that particles which, by chance, are far away from 
the center of the potential will experience a strong back- 
dragging force and will hence on average move further in 
the near future. As discussed in [|0| this effect can be 
also observed for the Rouse model of polymer dynamics. 
Note that on a qualitative level this type of heterogeneity 
is distinct from the type expressed by the roi-dependence 
of cTji j^(roi) for which the different mobilities are related 
to different environments. 

We gratefully acknowledge helpful discussions with S. 
Biichner and H.W. Spiess. 
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